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SUMMARY 


The technique of body-fitted coordinate systems, whereby curvilinear coor- 
dinate systems having coordinate lines coincident with all boundaries are gen- 
erated numerically as solutions of an elliptic partial differential system, is 
applied in numerical solutions of the complete time-dependent compressible and 
incompressible Navier-Stokes equations for laminar flow and to the time-depend- 
ent mean turbulent equations closed by modified Kolmogorov hypotheses for tur- 
bulent flow. Coordinate lines are automatically concentrated near to the 
bodies at higher Reynolds number so that accurate resolution of the large gra- 
dients near the solid boundaries is achieved. Two-dimensional bodies of arbi- 
trary shapes are treated, the body contour (s) being simply input to the pro- 
gram. The complication of the body shape is thus removed from the problem. 


INTRODUCTION 


The use of numerically generated boundary- fitted coordinate systems has 
made possible the development of numerical solutions of the Navier-Stokes equa- 
tions that can treat bodies of arbitrary shapes as easily as simple bodies. 

Codes can be written that are independent of the body or boundary shape, which 
may even be changing with time. 

These solutions are based on a method of automatic numerical generation of 
a general curvilinear coordinate system with coordinate lines coincident with 
all boundaries of a general multi-connected region containing any number of 
arbitrary shaped bodies (ref. 1). The curvilinear coordinates are generated as 
the solution of two elliptic partial differential equations with Dirichlet 
boundary conditions, one coordinate being specified to be constant on each of 
the boundaries, and a distribution of the other being specified along the bound- 
aries. Regardless of the shape, number, or movement of the bodies and regard- 
less of spacing of the curvilinear coordinate lines, all numerical computa- 
tions, both to generate the coordinate system and to subsequently solve the 
Navier-Stokes equations or any other partial differential equations on the 
coordinate system, are done on a rectangular grid with a square mesh, (i.e., in 
the transformed plane). The physical coordinate system has been, in effect, 
eliminated from the problem, at the expense of adding two elliptic equations to 
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the original system. Since the curvilinear coordinate system has coordinate 
lines coincident with the surface contours of all bodies present, all boundary 
conditions may be expressed at grid points, and normal derivatives on the bodies 
may be represented using only finite differences between grid points on coor- 
dinate lines, without need of any interpolation even though the coordinate sys- 
tem is not orthogonal at the boundary. Numerical solutions for the lifting and 
non-lifting potential flow about Karman-Tref f tz airfoils using this coordinate 
system generation show excellent agreement with the analytic solutions (refs. 1 
and 2). 

This method of automatic body-fitted curvilinear coordinate generation has 
been used to construct finite-difference solutions of the full time-dependent 
Navier-Stokes solutions for unsteady viscous flow about arbitrary two-dimension- 
al airfoils (refs. 2,3,4, and 5) and submerged and partially submerged hydro- 
foils (refs. 4 and 5). 

A method of controlling the spacing of the coordinate lines encircling the 
body has been developed in order to treat higher Reynolds number flow, since 
the coordinate lines must concentrate near the surface to a greater degree as 
the Reynolds number increases. The solution shows excellent comparison with 
the Blasius boundary layer solution for the flow past a semi-infinite plate 
(refs. 6 and 7). Solutions are also being developed for compressible viscous 
transonic flow with both subsonic and supersonic free streams, and for compres- 
sible turbulent flow. 


BOUNDARY-FITTED COORDINATE SYSTEMS 


The basic approach of constructing body-fitted curvilinear coordinate sys- 
tems in general multi-connected regions as the solution of an elliptic boundary 
value problem has been discussed in previous publications (refs. 1,3 and 8), 
and reference to related work by others has been made therein. A detailed re- 
port of the technique and code is now available (ref. 8), together with the 
code, a user's manual, and instructions with illustrations of its application to 
the numerical solution of partial differential equations. 

Certain considerations must be taken into account in the choice of a suit- 
able elliptic generating system for the coordinates as discussed in references 
3 and 8. The system chosen in the present work allows considerable control of 
the coordinate line spacing as is illustrated in reference 8. Control of the 
spacing of the coordinate lines on the body is easily accomplished, since the 
points on the body are input to the program. The spacing of the coordinate 
lines in the field, however, must be controlled by varying the elliptic gener- 
ating system for the coordinates. One method of variation is to add inhomogen- 
eous terms to the right sides of the Laplace equation, so that the generating 
system becomes 

5 + ? = P(5,n) , n + n = Q(5,n) (1) 

Since it is desired to perform all numerical computations in the uniform 
rectangular transformed plane, the dependent and independent variables must be 
interchanged in eq. (1). This results in the coupled system (see ref. 8 for 
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the transformation relations) . ^ 

axp- - 23x + yx = - J“^[x P(?,n) + x Q(5,n)l 

99 9H nri ? n 

+ yyj^f^ = - J^[y^P(5,n) + y^Q(5,n)] 

where 

2 2 2 2 
a = X + y 3 = x_x + y^y Y = + 7^ J = x^y - x y^ 

The inhomogeneous functions P(5,n) and Q(C>n) allow coordinate lines to be 
attracted to specified lines and/or points in the field or on the boundaries as 
discussed in detail in reference 8. It is necessary to give some consideration 
to the rapidity with which the spacing varies else truncation error effects in 
the form of numerical diffusion, possibly negative, may be introduced. All 
derivatives in equation (2) are approximated by second-order central finite 
difference expressions. The resulting difference equations are given in ref- 
erence 8. The set of nonlinear simultaneous difference equations is solved by 
point SOR iteration. 

In the present application to the Navier-Stokes equations, all diffussive 
space derivatives in the transformed equations are represented by second-order, 
central difference expressions. Both second-C"der central and. second-order 
upwind differences have been investigated for convective derivatives. Deriva- 
tives off a boundary are represented by backward difference expressions being 
solved simultaneously by point SOR iteration at each time step. 


(2a) 

(2b) 


LAMINAR FLOW ABOUT MULTIPLE AIRFOILS 


With the velocity and pressure as the dependent variables the transformed 
Navier-Stokes equations are (ref. 8) 

Uj. + [y^(u^)^ - y^(u^)^]/J + tXg(uv)^ - x^(uv)^]/J 


+ (Qu^ + Pu^)/R 

V + [y (uv) - yp(uv) ]/J + [x (v^) - x (v^) ]/J 

t ri99n 9nri9 

+ (x_p - x p_)/J= (av^_ - 23v_ + yv )/RJ, 


+ (Qv^ + Pv^)/R 

■ ■*'nn - 

- - x^u^Ky^v^ - y^v^) - (x^v^ - x^v^) - 


ap^P - 23p^^ + YP„„ + (Qp„ + ppf)J2 = - (y^^u^ - y^V"" 


(3) 


(4) 


(5) 
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where R is the Reynolds number 


(y + x^v - X v^)/J (6) 

Equation (5) is the transformed Poisson equation for the pressure, obtained by 
taking the divergence of the Navier-Stokes equations. 

The boundary conditions are 

u = V = 0 on body surface (7a) 

u = COS0, V = sin0, p = 0 on remote boundary (7b) 

The pressure at each point on the body was adjusted at each iteration by an 
amount proportional to the velocity divergence evaluated using second-order one- 
sided differences for the n-derivative on the body. 

The body force components are obtained from 
and shear forces around the body surface: 

2 

Fjj = + py^ ^ dC 

2 

Fy = - / Px^ d? - - ^ d? 

and the lift and drag coefficients are given by 

C_ = F cos0 - F sin0 
L y X 

C = F sin0 + F COS0 
D y X 

where 0 is the angle of attack. 

In the velocity-pressure formulation it is necessary to calculate the body vor- 
ticity before applying equation (8) from 

Figure la shows the coordinate system for a multiple airfoil consisting of two 
Karman-Tref f tz airfoils, one simulating a separated flap. Coordinate system con- 
trol was used to attract the coordinate lines strongly to the first ten lines 
around the bodies and to the intersections of the cut between the bodies with 
the trailing edge of the forebody and the leading edge of the aftbody. Veloc- 
ity vectors for the viscous flow at R - 1000 are shown in figures lb - e. In- 
cluded in these plots are detail views of the slot between the airfoils, the 
trailing edge of the aft airfoil, and the separated region about the aft airfoil. 

Results have recently been obtained by Hodge at Flight Dynamics Laboratory, 
Wright-Patterson Airforce Base, (ref. 5), for a single NACA 0018 airfoil at R = 
41,400 which compare very well with experimental drag values. 


the integration of the pressure 

(8a) 

(8b) 

(9a) 

(9b) 


TURBULENT SHEAR FLOW AROUND A CIRCULAR CYLINDER 


The aim of the following analysis Is to select a turbulence model for the 
computation of mean turbulent flow fields around two-dimensional bodies of arbl- 
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trary shapes in the body-fitted coordinate systems. As a test case the coordi- 
nate system for a circular cylinder was generated over which the model equations 
of mean turbulent flow were solved. 

As is. well known, the mean turbulent flow equations and the moment equa- 
tions of the desired order can be obtained from the non-steady three-dimension- 
al Navier-Stokes equations. The specification of the unknown correlations forms 
the closure problem. This line of approach for closing the system of equations 
was initiated by Chou (ref. 9) and Rotta (ref. 10). Since then much work has 
been done on the closure problem (see refs. 11 and 12). 

A philosophically different approach was initiated by Kolmogorov (ref. 13) 
in which the turbulent transport equations were written down heuristically to 
model the physics rather than to model the unknown correlations resulting from 
the averaged Navier-Stokes equations. A set of equations similar to those of 
Kolmogorov were proposed by Saffman (refs. 14 and 15). The prediction capabil- 
ities of the Kolmogorov’s equations have not been investigated at all, while 
Saffman himself has used his equations only in the cases of plane Couette flow, 
plane Poiseuille flow, and the two-dimensional jet and wake. 

The- turbulence model chosen in this paper to describe the two-dimensional 
mean turbulent flow around finite bodies is comprised of the energy transport 
equation of Kolmogorov’s model and the vorticity-density transport equation of 
the Saffman model. Both models have been described by Saffman O^ef. 15), 


Using Cartesian tensor rotation and the summation convention the equations 
of continuity and momentum for an incompressible flow are 


BU. 
. 1 

Bt. 


BU. 

1 

Bx . 
3 

BU. 


= 0 


+ u 


j 3x 


- + vvju. - ~ 

p 3x. 1 1 3x. 

1 3 




( 11 ) 


( 12 ) 


where t^ is the physical time, p .the density, v the kinematic viscosity, U. the 
Cartesian components of the mean velocity vector and u. the fluctuating compo- 
nents. An overbar denotes the average. 


To close the system of equations (11) and (12), the eddy viscosity formu- 
lation due to Boussinesq is used 


2 -c 


where 


u.u. = *— e6. . - 2v„S. . 
1 j 3 xj T JLj 


6. . is the Kronecker delta 
. 3U. 3U. 

3 1 


(13) 


e = I u.u^ 


the mean turbvilence energy 
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= kinematic eddy viscosity 

The Kolmogorov-Saffman model now depends in specifying as 

a^e 

Vt = -f- (14) 

where a is a constant and 0 is the square root of the turbulence vorticity- 
density defined as 


1/2 


1/2 


0 = (<l>) = 

oj. being the vorticity fluctuations. The quantities e and <j) are to be deter- 
mined from their own transport equations which are 


^ + V - + if: 

1 J J 3 


(15) 


1 J J J 

where — 1/2 , 

P = -;r I \ K — J J » “ » h constants 

2 dX. 0 


(16) 


Q = 2(s. )■ 

2 32 


1 3x.0x, 

3 J 

The boundary conditions for equations (11), (12), (15), and (16) are 
at the wall: S^u.^ 

u. = 0, 1 - 0, ? . ' ° . .'y 


(17) 


at infinity: 


U. = U., e = oore, <b = 0 

X “X oo ^ 


where u. is the friction velocity and S a non-dimensional constant having the 

value o? about 8000 and a =0.3, b=5/3. 

o 

The required equations to predict the mean turbulent flow in two dimensions 
are equations (11)- (16). For the numerical computation the vorticity-s tream 
function formulation i^ adopted in which 

^1 3X2 ’ ^ ’ 

SO that in place of (11) and (12) the equations are 


and 


= - 


. TT „ 2— . 

-r— + U . -r = vV W + ^ 

3 t^ J 3 x . 1 3 x, 3 x 


(Ui -U 2 ) + ( — 


'1"^2 


3x, 


3x. 


-j) 


(18) 

(19) 
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The pertinent equations are now non-dimensionalized as follows 

E = e/V^, 0) = Lw/V, X = L^?/V^K^, 

T = V^/2VL, i|; = '5i/VL, t = Vt^/L, P = LP/V, 

Q = L^Q/V^, R = VL/v, X = x^/L, y = 

2 2 

u = Uj^/V. V - Uj/V, 72 - 

where V and L are the characteristic velocity and length respectively. 


(20) 


In the expression for x the function K is assumed to be a non-dimensional 
function of the friction velocity u.. This function has been introduced to absorb 
the variation of (p at the surface as indicated in equation (17) so as to make 
X = a constant at the wall. Now at the wall 


X == a 

u. ^ o 

where U. = so that the choice x = 1 at the wall gives 

V O ^ N 

K = a SR (21) 

o *5? 

From the analysis of Saffman (ref. 15) the value of S for smooth walls is about 
8000. I Since R is the free stream Reynolds number the numerical values of K are 
quite large. Neglecting the gradients of K in the x equation, the non-dimension- 
al equations for numerical computation are 


= -(jO 

X + - t \ =-~v2x+ M(T X -H T X )+AX + B 

t y X X y H XX y y 


( 22 ) 

(23) 


where 


H = 


T = 


R 

l+CR T 

o E 

1/2 


2Kx 

and the variable subscripts denote partial differentiations. The 
B, M and C depend on the choice of the surrogate variable X which 
below. 

a = (jj: a=zv“i.D = Hi __ 

XX y yy X xy y 

2. 

B = 2TQ, M = 1, C = 1 


oj: 

A = 2v2t, 


M = 4, C ^ 


* ^ 1/2 

E: 

A =-Kx 

X5 

A = P-bKx 


1/2 


(24) 


values of A, 
are described 


(25) 

(26) 

(27) 


The boundary conditions are 


at the body surface: ip=0, X“l> E=0, 

at infinity: oH-0, X"^> E->0. (28) 

Equations (22) and (23) have been transformed to the body-fitted coordinates and 
subjected to the finite difference approximations similar to that used in the 
laminar calculation. Velocity vector plots and the energy distribution for a cir- 
cular cylinder at R = 5x10^ are shown in figures 2a and 2b. In figure 2b 


Jlnr 


.0576n(1.2)'^ 


r = non-dimensional radial distance. 
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Figure 2. Turbulent flow past a circular cylinder (R = 500,000, 
t = 0.29; transition introduced at t = 0. 12) :(a)velocity vectors 
from 90® to the aft stagnation point, (b) turbulence energy 
distribution (l. front stagnation, 2. aft stagnation, 3. 90®) 


